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SUMMARY.  The  numerical  results  of  a  class  of  problems  of  linear 
elastic  stability  problems  subjected  to  nonconservative  forces  and 
under  various  support  conditions  are  presented  here.  A  single  solu¬ 
tion  formulation  by  which  these  results  have  been  obtained  is  descri¬ 
bed.  Accuracy  of  these  results  compared  with  those  reported  in  the 
literature  is  discussed. 

I.  INTRODUCTION.  Any  particular  subject  of  investigation  in 
applied  sciences  is  always  motivated  by  the  desire  to  understand  some 
natural  phenomena  and  hopefully  to  utilize  the  results  of  such  an 
investigation  for  the  benefit  of  human  activities.  The  study  of 
structural  behavior  under  nonconservative  loads  is  of  no  exception. 

Since  follower  forces  are  a  special  class  of  nonconservative  forces  [1] , 
one  is  surprised  to  encounter  frequently  the  question  as  to  the  relation 
between  such  a  study  and  a  real  engineering  problem.  Physically,  a 
follower  force  is  simply  one  whose  direction  follows  the  structural 
deformation  as  in  comparison  with  a  dead  load  which  acts  in  a  fixed 
direction  independent  of  deformation.  Some  obvious  examples  of 
follower  forces  are:  thrust  at  the  tail  of  a  flexible  rocket,  jet 
engine  thrust  of  an  airplane,  thrust  on  the  propeller  shaft  of  a  ship, 
etc.  Other  examples  such  as  the  pressure-and-curvature  induced  forces 
included  in  the  gun  dynamic  studies  are  less  obvious  [2] 

Since  the  problems  of  follower  forces  are  non-self-adjoint  their 
treatment  is  more  difficult  than  that  for  the  self-adjoint  problems. 

In  the  classical  paper  by  Beck  [3],  it  was  demonstrated  that  the 
stability  nature  of  a  nonconservative  problem  can  be  quite  different 
than  that  of  a  conservative  one.  For  these  reasons,  a  systematic 
approach  to  this  class  of  problems  and  an  understanding  of  some  of  the 
basic  problems  involving  follower  forces  are  desirable. 


The  purpose  of  this  paper  is  to  present  a  single  solution  approach 
to  a  class  of  problems  of  follower  forces,  including  several  classical 
examples,  to  present  the  numerical  results  so  obtained,  and  to  discuss 
the  accuracy  compared  with  those  already  published  in  literature. 
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In  Section  II,  the  class  of  problems  will  be  defined  by  a  general 
form  of  a  differential  equation  and  a  set  of  boundary  conditions.  The 
solution  formulation  and  its  basis  aregiven  in  Section  III.  Numerical 
results  of  some  specific  problems  are  given  in  Section  IV  together  with 
a  discussion  and  comparisons  with  data  available  in  literature. 

II.  A  CLASS  OF  PROBLEMS  SUBJECTED  TO  FOLLOWER  LOADS.  The  class  of 
problems  considered  in  this  paper  can  be  described  by  the  differential 
equation 


y""  +  P(x)y"  +  X2y  =  0  (1) 

where  y(x)  denotes  the  lateral  disturbance  of  a  beam,  as  a  function  of 
the  abscissa  x,  P(x)  is  the  axial  force  always  tangent  to  the  deformed 
axis,  and  X  is  the  eigenvalue.  As  usual,  a  prime  denotes  differentiation 
with  respect  to  x. 

Eq.  (1)  is  a  non-self-adjoint  differential  equation  (thus  noncon¬ 
servative  problem)  except  for  P(x)  =  constant.  If  the  axial  force  P(x) 
remains  fixed  in  the  direction  of  the  undeformed  axis,  the  problem  would 
be  of  conservative  nature  and  the  differential  equation  a  self-adjoint 
one. 


y""  +  .[P(x)y']'  +  A2y  =  0  (1') 

Both  Eqs .  (1)  and  (1')  are  well  known  and  the  derivations  are  simple  and 
they  follow  the  procedures  given  in  such  textbooks  as  that  by  Timoshenko 
and  Gere  [4].  Boundary  conditions  considered  will  be  in  the  following 
form: 


y"«  (0)  +  P (O)y'  (0)  +  k1(0)  =  0  (2a) 

-y" (0)  +  k2y'(0)  =  0  (2b) 

-  y"'  (1)  -  (l-k5)P(l)y'(l)  +  k3y(l)  =  0  (2c) 

y" (1)  +  k4y' (1)  =  0  (2d) 

where  kj_,  k2  are  the  deflection  and  rotation  spring  constants  at  x  =  0 

and  k3,  k4  are  the  same  at  x  =  1.  The  constant  k5  is  related  to  a 
"constant  of  tangency"  Kq  by  equation 

K0  =  k5  _  1  ^ 
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so  that  Eq.  (2c)  becomes 

-y'"  (1)  +  K0P(l)y'(l)  +  k3y(l)  =  0  (2c ' ) 

where  now,  if  P(l)  /  0,  0  =  Kgy ' ( 1)  denotes  the  angle  that  P(l)  is  to  be 
rotated  with  respect  to  the  tangent  of  the  beam  at  x  =  1  (Figure  1). 

Eqs.  (2)  simply  state  that  the  total  shear  force  and  moment  at 
x  =  0  and  x  =  1  must  be  zero.  As  kj  approaches  infinity,  Eq.  (2a) 
requires  that  y(0)  =  0.  Thus  a  zero  deflection  boundary  condition  is 
arrived  at.  Similar  options  are  provided  for  by  other  spring  constants 
k2,  kj  and  k4. 

Three  different  P(x)  will  be  considered  in  this  paper:  (1)  P(x)  = 

P,  a  constant,  (2)  P(x)  =  q(l-x),  and  (3)  P(x)  =  qQ/2(l-x)2  where  P 
represents  a  concentrated  force  at  x  =  0,  q  is  a  uniformly  distributed 
follower  force  density  and  q0  denotes  the  maximum  of  a  linearly  varied 
follower  force  density.  With  the  special  boundary  conditions  of  a 
cantilever,  case  (1),  (2),  and  (3)  become  the  classical  problems  first 
solved  by  Beck  [3],  Leipholz  [5],  and  Hanger  [6],  respectively. 

III.  SOLUTION  FORMULATIONS.  The  solution  method  used  here  is  the 
finite  element  unconstrained  variational  formulation  which  has  proved 
to  be  efficient  and  simple  to.  use  for  solutions  of  non-self-adjoint 
problems  [7,8].  Finite  elements  are  used  in  the  usual  sense  that  the 
unknown  function  is  approximated  by  piecewise  cubic  splines.  An  uncon¬ 
strained  variational  statement  is  established  and  used  so  that  none  of 
the  boundary  conditions  need  to  be  satisfied  a  priori.  An  outline  of 
the  formulation  will  be  given  here. 

Introducing  an  adjoint  field  variable  y*(x),  it  is  a  simple  matter 
to  see  that  the  following  variational  statement  will  lead  to  the  differ¬ 
ential  equation  (1)  and  boundary  conditions  (2) : 

6I(y>y*)  =  0  (4a) 

1 

I  =  /  (y"y*"-P(x)y,y*'-P, (x)y*y*+X2yy*)dx 
0 

+  kiy(0)y*(0)  +  k2y'(0)y*'(0)  +  k3y(l)y*(l)  +  k4y ' (l)y* ' (1) 

+  k5P(l)y'(l)y*(l)  .  (4b) 

The  fact  that  Eqs.  (4)  lead  to  the  given  differential  equation  and 
boundary  conditions  for  y(x)  independent  of  y*(x)  implies  that  one  can 
take  the  variation  of  I  at  y*(x)  =  0  and  (<Sl)y*=g  =  0  still  leads  to 
the  original  problem.  Hence  our  formulation  begins  with 

(SI)y,E0  -  0  (5a> 
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or,  1 

/  [y"5y*"  -  P (x)y'6y*1  -  P(x)y'6y*  +  X2y<$y*]dx 
0 

+  k1y(0)6y*(0)  +  k2y '  (0)  <Sy* '  (0)  +  k3y(l)Sy*(l)  +  k4y'  (l)<5y*'  (1) 

+  k^y'  (1)  6y*  (1)  =  0  (5b) 

Finite  element  discretization  enters  when  the  beam  is  divided  into 
L  equal  elements  and  Eq.  (5b)  is  written  as 


l  J  [L3y('1')"<5y*('1-)"  -  LP(l)  (£)y^  '<Sy*(l) ' 
i=l  0 

-  LpW'y(i)V(i)  +  Aly(i)6y*(i)]d^ 

L 


+  kp^1)  (0)<5y*(1)  (0)  +  k2L2y^1^  *  (0)  6y*  ^ '  (0) 

+  k3y(L^>  (l)6y*(L)  (1)  +  k4L2y  '  (1)  6y*  '  (1) 

+  k5y(L^' (l)6y*(L) (1)  =  0  (6) 

In  obtaining  Eq„  (6)  from  (5b),  one  has  effected  a  change  of  coordinates 
from  x  (global)  to  £  (local)  such  that 

?  =  =  Lx-i+1 

d£  =  Ldx 

(i)  ^ 

y(x)  =  yUJ(S) 

y'  ^  =  ~h  y('x-)  =  L  ^  y('1’)  ^  =  Ly('1’) '  ^ 

etc. 


Introducing  generalized  coordinates  vector  and  shape  function  vector 

a (5)  such  that 


y(i,«)  -  aT(Ov(1) 

(8a) 

with 

yti)T  „  {Yi(i)  y2Ci)  y3Ci)  y^Dj 

(8b) 
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a(5)  = 


0 

-3 

2~ 

[f1  \ 

0 

1 

-2 

1 

u 

0 

0 

3 

-2 

5* 

0 

0 

-1 

1 

[s3, 

(8c) 


where  a  superscript  T  denotes  the  transpose  of  a  matrix-  One  observes  that 
Y1(l)  =  y^  (0)  ,  Y2(i)  =  y(i)'(0) 

Y3(i)  =  y(i)(l)  ,  Y4W-y(i)'(l)  (8d) 

The  counterparts  for  y*(£)  can  be  similarly  defined; 

In  terms  of  Y^,  Y*^,  a.,  Eq.  (6)  can  be  written  as: 

l  SY*WT{L3J  a"(Oa"T(S)dS  -  L /  P(i)  (S)a'  (£)a'T(Qd5 
i=0  0~  0  ~ 

-  L /  P(l),(5)a(C)a»T(?)d^  +  /  a(OaT(5)dC}Y(i) 

0  ~  L  0"  ~ 

+  6Y*<-1-)T{k.a(0)aT(0)  +  k„L2a'  (0)a,T(0)  }Y^ 

+  6Y*(L)T{k3a(l)aT(l)  +  k4L2a' (l)a,T(l)  +  k5a(l)a'T(l) }y(L)  *  0  (9) 

It  will  be  convenient  to  define  the  following  matrices: 


1  t  1  T 

A,  =  /  aCOa^Od?  ,  A2  =  /  a'COa'^Odf: 

0  0 

A  =  /  a»CC)a"T(Od?  ,  A  -  /  a(£)a'T(S)dS 
~3  0  ~4  0 
1  T  1 
A5  =  /  £a'CC)a'T(OdS  ,  A&  =  /  a(^a'T(?)dC 

0  0 

1  T 

A7  =  /  £V(5)a'  (Qd? 

0 

Bx  =  a(0)aT(0)  ,  B2  =  a'(0)a'(0) 

B3  =  a(l)aT(l)  ,  B4  =  a'(l)a'T(l) 


(10) 


B5  =  a(l)a'T(l) 
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In  terms  of  the  matrices  defined  in  (10) ,  Eq.  (9)  is  written  as: 

L  2 

V  6Y*^1^)T{L3A_  +  —  A.  -  LM  JY^ 

. u ^  ~3  T  ~1  ~p  ~ 

i=l  L  r 

+  6Y*(1)T{k1B1  +  k2L2B2}Y(1) 

+  SY*(L)T{k3B3  +  k4L2B4  +  k5B5}Y(L)  =  0 


(11) 


where  the  matrix  M_  is  defined  as 

t  <  i  i 

M  =  /  P(l)a'a'TdC  +  /  P(l)'aa'Td?  .  (12) 

~P  0  ~  ~  o 

To  proceed  further,  it  is  necessary  to  know  the  specific  form  of  P(x). 

As  we  have  mentioned  earlier,  three  different  forms  of  P(x)  will  be 
considered. 


CASE  I.  P(x)  =  P,  a  Constant.  In  this  case,  one  has 

P(x)  -  P(l)(0  =  P 


P'(x)  =  LP(l),(£)  =  0 


Thus 


M  =  P /  a'a'd E,  =  PA 


"2 


(13) 


(14) 


CASE  II.  P(x)  =  q(l-x). 

P(x)  =  P(l)(0  =  ^  (L-i+1-0 

P^'  (O  =  -  1 
L 


Thus, 


„  1  T  1  T 

M*  f  {[L  -  (i— 1) ] /  a'a»*d5  -  /  Ca'a^dC 

~P  "  0~~  0~~ 


or 


Mp  =a([L  -  (i-l))A2  -  a5} 


(15) 


(16) 
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CASE  III.  P(x)  =  q0/2(l-x)2o 


Thus , 


or. 


P(i)(5)  =  (L-i+1)  2  -  2 (L-i+1) E,  +  C2] 

Z  L) 

^  [(L-x+i)  -  q 


M_  =  {(L-i+1)2/  a'a'TdC  -  2(L-i+l)/  ?a'a'Td£ 

-P  2LZ  o  -  ~  0  ^  ~ 

1  T 

+  /  ?2a'a'  d£) 

0  ~  ~ 

1  1 

-  3°.  {(L-i+1)/  a'a,Td£  -  /  5a'a,Td£;} 

L2  0“  ~  0  ~  ~ 


{  (L-i+l)2A2 


2 (L-i+1) Ag  +  A?) 


(17) 


•  .  JL  { (L-i+1)  A  -  A  }  .  (18) 

L2  ~2 

With  M_  defined  for  all  three  cases  in  Eqs.  (14),  (16),  and  (18)  respec¬ 
tively,  one  can  now  assemble  Eq.  (11)  into  a  global  matrix  equation. 
Introducing  the  global  generalized  coordinate  vectors  Y  and  Y*  as: 

yT-tytO  y,(»  y,™  y<»  v'2>  y<2> . y'L’  r  'l>  ) 

-  I  Z  3  4  5  4  4  Q9) 

Y*T  =  {y  Y  Y  Y  *  Y  *  Y  *  (2) . Y  Y  *  (k) } 

~  123434  34 

Eq.  (11)  now  can  be  written  in  terms  of  Y  and  6Y*  as 

6Y*T{K  -  \2M}Y  =  0  (20) 

a/ 

where  the  global  matrices  K  and  M  are  formed  by  properly  placing  the 
local  matrices  defined  in  Eqs.  (10)  according  to  the  correspondence 
between  the  local  and  global  generalized  coordinates  indicated  in  Eqs. 
(19).  Now  since  6Y*  are  not  subject  to  any  constraint  conditions,  Eq. 

(18)  reduces  to 


(K  -  X2M)Y  =  0 

which  is  solved  for  the  eigenvalue  X  and  the  eigenvector  Y. 


(21) 
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IV.  NUMERICAL  RESULTS  AND  DISCUSSION.  It  is  yell  known  that  the 
eigenvalue  X  dictates  the  stability  of  the  column:  a  pure  imaginary 
X  is- associated  with  a  stable  vibration,  a  real  X  with  instability  of 
divergence,  and  a  complex  X  with  instability  of  flutter  [10]. 


Only  cantilevered  columns  will  be  considered  here.  It  will  be  seen 
that  in  all  three  loading  cases,  the  cantilevered  columns  reaches  an 
instability  condition  of  flutter„ 

CASE  I.,  P(x)  =  P  =  Constant.  The  characteristic  equation  in  close 
form  was  obtained  by  Beck  [3]  as 

2X2  +  Q2  +  2X2coshacosB  +  QXsinhasinB  =  0  (22) 

where 

a2  =^X2  +  £  +  8- 

4  2 

_  (23) 

.  6*  jfcTf  ♦  I 

For  a  given  Q,  the  eigenvalue  X  can  be  calculated  from  Eq.  (22)  and  there 
are  an  infinite  number  of  X  solutions  for  each  Q.  Eq.  (22)  is  solved  for 
two. lowest  branches  of  X  using  an  iterative  procedure.  The  results  are 

given  in  Table.  I.  The  critical  load  thus  obtained  is 

,  '•  | 

Qcr  =  2.0318TT2  =  20.053 

which  agrees  well  with  the  value  obtained  originally  by  Beck  a?  QqR  = 
20.05.  The  results  for  four  lowest  eigenvalues  presently  obtained  using 
our  finite  element-unconstrained  variational  formulations  are  also  shown 
in  Table  I.  The  first  two  branches  obviously  agree  well  with  those  from 
the  exact  characteristic  equation.  It  should  be  pointed  out  that  the 
numerical  solutions  to  the  Beck  problem  given  in  Reference  [4]  appear 
to  be  inaccurate.  A  plot  of  the  eigenvalue  curve  showing  the  coalescense 
of  the  two  lowest  branches  is  given  in  Figure  2.  The  data  from  Reference 
[4]  are  indicated  by  small  circles »  The  fact  that  these  data  points  do 
not  fall  on  a  smooth  curve  further  add  to  the  doubt  about  their  accuracy. 
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CASE  II.  P(x)  =  g(l-x) .  The  numerical  values:  of  the  four  lowest 
eigenvalues  up  to  the  first  critical  load  are  giyen  in  Table  II  which 
the  critical  load  is  shown  to  be 

qCR  =  4.05911T2 

compared  with  data  given  by  Leipholz  as  4.1238tt2  =  40.7  [11]  and  again 
as  4.2058'rr2  =  41.51  [12].  The  coalescence  of  the  first  two  branches 
of  eigenvalues  is  again  shown  in  Figure  2. 

CASE  III.  P(x)  =  q0/2(l-x)2.  Similar  data  for  this  case  are 

presented  in  Table  III  and  in  Figure  3.  The  critical  load  of  flutter 
is  obtained  as 

qoCR  =  15. 2687ir2 

In  comparison,  the  value  obtained  by  Hauger  was  qQcR  *  158.2  -  16.09^tt2 
[6]  and  that  by  Leipholz,  qoCR  =  150.80  =  15.279tt*  [12]. 
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TABLE  I.  NUMERICAL  VALUES  OF  FIRST  TWO  LOWEST  EIGENVALUES  OF 
A  CANTILEVERED  COLUMN  WITH  P(x)  =  P,  A  CONSTANT 
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TABLE  II.  NUMERICAL  VALUES  OF  FIRST  FOUR  LOWEST  EIGENVALUES  OF 
_ A  CANTILEVERED  COLUMN  WITH  P(x)  =  q(l-x) _ 
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TABLE  III.  NUMERICAL  VALUES  OF  FIRST  FOUR  LOWEST  EIGENVALUES  OF 
A  CANTILEVERED  COLUMN  WITH  P(x)  =  qQ(l-x)2/2 


V© 

© 

V' 

© 

00 

r^ 

*fr 

LO 

© 

00 

00 

00 

CM 

Tt 

© 

iH 

• 

• 

• 

• 

• 

LO 

rH 

f-H 

CM 

iH 

rH 

rH 

rH 

LO 

rH 

rH 

O 

00 

t^ 

r- 

to 

• 

CM 

CM 

rH 

o 

rH 

rH 

CM 

rH 

• 

• 

• 

• 

v© 

to 

rH 

LO 

rH 

rH 

KV 

o 

CM 

o 

H 

00 

rH 

• 

00 

LO 

LO 

oo 

LO 

<N 

O 

• 

• 

• 

• 

00 

r- 

VD 

rH 

.  LO 

rH 

rH 

o 

© 

L0 

© 

o 

L0 

CM 

LO 

• 

rH 

rf 

LO 

O 

to 

LO 

© 

• 

• 

• 

t 

O 

Gi 

00 

CM 

L0 

rH 

rH 

o 

v© 

cn 

© 

, 

© 

LO 

o 

o 

rH 

to 

CM 

to 

o 

r^ 

o 

• 

# 

« 

• 

to 

CM 

rH 

rH 

cm 

SO 

CM 

rH 

N 

N. 

rH 

r< 

to 

Tt 

r< 

<y 

13 


Figure  1.  Boundary  Condition  Associated  with  a  Follower  Force: 
Constant  of  Tangency  Kq„ 
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Figure  2„  Two  Lowest  Eigenvalues  vs.  Load  Parameters  Q/tt2  and  q/ir2. 
A  Cantilevered  Column  of  Case  I  and  II. 

(Data  shown  in  dots  (®)  are  from  Reference  [4]) 
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Figure  3»  Two  Lowest  Eigenvalues  vs.  Load  Parameter  q0/Tt 
A  Cantilevered  Column  of  Case  III. 
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